### SCRIPT FOR:
###
### World Population Growth Over Millennia: 
### Ancient and Present Phases with a Temporary Halt In-Between
###
### VERSION: 2023-JAN-26


rm(list=ls())
gc()

setwd("C:/Users/miros/Desktop/Research/World Population Growth")



# LIBRARIES ---------------------------------------------------------------
library(readxl)
library(scales)
library(ggplot2)
library(ggpubr)
library(grid)




# DATA --------------------------------------------------------------------

# __ HYDE 3.2: History Database of the Global Environment -----------------
#    NOTE: The data for 10 000 BCE to 2015 CE are taken from HYDE 3.2. Data 
#          were downloaded from https://doi.org/10.17026/dans-znk-cfy3 
#          [accessed 11-MAY-2022]. This research uses population 'TOTALS'
#          included in files 'popc_r.txt' available separately in the 
#          folders 'baseline', 'lower', and 'upper'.
#          For details about the dataset construction, see associated article: 
#          Klein Goldewijk, K., Beusen, A., Doelman, J., and Stehfest, E. 
#          (2017): Anthropogenic land use estimates for the Holocene - 
#          HYDE 3.2, Earth Syst. Sci. Data, 9, 927-953 
#          at https://doi.org/10.5194/essd-9-927-2017.

hyde <- read_excel("01-data/HYDE 3.2/HYDE 3.2_popc_r_aggregate.xlsx", 
                   sheet = "combined")

# Adding a note on source
hyde$source <- c("HYDE 3.2")





# __ Gapminder (version 6.0) ----------------------------------------------
#    NOTE: The data for 1800 to 2100 are taken from Gapminder 6.0. Data were
#          downloaded from online spreadsheet at https://gapm.io/d_popv6
#          [accessed 11-MAY-2022].
#          For details on dataset construction see: 
#          https://www.gapminder.org/data/documentation/gd003/

gapminder <- read_excel("01-data/Gapminder 6.0/GM-Population - Dataset - v6.xlsx", 
                        sheet = "data-for-world-by-year")

# Erasing irrelevant variables
gapminder$geo <- NULL
gapminder$name <- NULL

# Harmonizing variable names
names(gapminder)[names(gapminder) == "time"] <- "year"
names(gapminder)[names(gapminder) == "Population"] <- "baseline"

# Adding empty columns (to make the object consistent with other data objects)
gapminder$lower <- NA
gapminder$upper <- NA

# Adding a note on source
gapminder$source <- c("Gapminder 6.0")





# __ UN World Population Prospects (2019) ---------------------------------
#    NOTE: The data from UN World Population Prospects (2019) consist of
#          two parts: (1) Population estimates for the 1950-2020 period
#          and (2) Population projects for the 2020-2100 period. Both parts
#          are loaded separately. Data were downloaded from the UN depository
#          at https://population.un.org/wpp/Download/Standard/Population/
#          [accessed 11-MAY-2022]. Details about the data can be found on
#          the same webpage.


# ____ UN 1950-2020 population estimates ----------------------------------
UN_1950_2020 <- read_excel("01-data/UN World Population Prospects 2019/WPP2019_POP_Totals_F01_1_TOTAL_POPULATION.xlsx", 
                            sheet = "1950 to 2020")

# Harmonizing variable names
names(UN_1950_2020)[names(UN_1950_2020) == "Type"] <- "year"
names(UN_1950_2020)[names(UN_1950_2020) == "World"] <- "baseline"

# Adding empty columns (to make the object consistent with other data objects)
UN_1950_2020$lower <- NA
UN_1950_2020$upper <- NA

# UN population prospects are in thousands, so they are multiplied by 1000
# to get the (full) values consistent with other data sources
UN_1950_2020$baseline <- UN_1950_2020$baseline*1000

# Adding a note on source
UN_1950_2020$source <- c("UN Population Prospects (2019)")




# ____ UN 2020-2100 population prospects ----------------------------------
UN_2020_2100 <- read_excel("01-data/UN World Population Prospects 2019/WPP2019_POP_Totals_F01_1_TOTAL_POPULATION.xlsx", 
                           sheet = "2020 to 2100")

# Harmonizing variable names
names(UN_2020_2100)[names(UN_2020_2100) == "Year"] <- "year"
names(UN_2020_2100)[names(UN_2020_2100) == "Medium variant"] <- "baseline"
names(UN_2020_2100)[names(UN_2020_2100) == "Low variant"] <- "lower"
names(UN_2020_2100)[names(UN_2020_2100) == "High variant"] <- "upper"

# UN population prospects are in thousands, so they are multiplied by 1000
# to get the (full) values consistent with other data sources
UN_2020_2100$baseline <- UN_2020_2100$baseline*1000
UN_2020_2100$lower <- UN_2020_2100$lower*1000
UN_2020_2100$upper <- UN_2020_2100$upper*1000

# Erasing year 2020 which is already included in the 1950-2020 estimates
UN_2020_2100 <- subset(UN_2020_2100, year != 2020)

# Adding a note on source
UN_2020_2100$source <- c("UN Population Prospects (2019)")





# __ Taagepera (1979) ----------------------------------------------
#    NOTE: The data collected in Taagepera, Rein. 1979. "People, 
#    Skills, and Resources: An Interaction Model for World Population 
#    Growth." Technological Forecasting and Social Change 13(1): 13-30.
#    https://doi.org/10.1016/0040-1625(79)90003-9. It relies on multiple
#    resources that are directly referred to as 'source'.

taagepera <- read_excel("01-data/Taagepera 1979/estimates_from_literature.xlsx")

# Replacing the source information with 'Taagepera (1979)'
taagepera$source <- "Taagepera (1979)"





# __ Miscellaneous sources ------------------------------------------------
# NOTE: This dataset is a collection of sources summarized on Wikipedia
#       https://en.wikipedia.org/wiki/Estimates_of_historical_world_population#cite_note-26
#       and then manually reshaped into the format loaded here.

misc <- read_excel("01-data/Wikipedia sources/pop_estimates_review.xlsx")





# __ MERGE of objects -----------------------------------------------------
# NOTE: The four data files are merged into a single data object which is 
#       used from now on.

population <- rbind(hyde,
                    gapminder,
                    UN_1950_2020,
                    UN_2020_2100,
                    taagepera,
                    misc)

# Getting rid of irrelevant objects
rm(hyde, gapminder, UN_1950_2020, UN_2020_2100, taagepera, misc)




# FUNCTION ----------------------------------------------------------------
# Defining a function that can reverse and log the scale
# NOTE: 'ggplot2' can handle only one type of scale transformation at the
#       same time. However, this visualization requires logging and
#       reversing the scale at the same time. Therefore, this function is
#       defined to accommodate both scale transformations. It is used later
#       in the script during the visualization. 

reverselog_trans <- function(base = exp(1)) {
  trans <- function(x) -log(x, base)
  inv <- function(x) base^(-x)
  trans_new(paste0("reverselog-", format(base)), trans, inv, 
            log_breaks(base = base), 
            domain = c(1e-100, Inf))
}





# ANALYSIS ----------------------------------------------------------------

# Subtracting 2200 from the year variable so it is possible to use logarithm
population$year_min2200 <- 2200-population$year

# Defining the functions that will be visualized
f.PT.new.1 <- function(x) (3.82*10^3)/log(1.25+exp((1980-(2200-x))/25.5))^0.716
f.PT.new.2 <- function(x) (3.82*10^3)/(((1980-(2200-x))/25.5))^0.716
f.PT.old.1 <- function(x) (2.3*10^3)/log(34000+exp((100-(2200-x))/25.5))
f.PT.old.2 <- function(x) (2.3*10^3)/((100-(2200-x))/25.5)




# __ FIGURE 1: Present period -----------------------------------------------

ggplot(subset(population, year >= 400), 
       aes(x = year_min2200, y = baseline/1000000)) +
  geom_point(size = 1) +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.005, size = 0.2) +
  scale_x_continuous(breaks = c(2200, 2100, 2000, 1900, 1800, 1700, 1600, 1500, 1400, 1300, 1200, 1100, 1000, 900, 800, 700, 600, 500, 400, 300, 200, 100),
                     labels = c("0(1)\n", "", "", "", "", "500\n", "", "", "", "", "1000\n", "", "", "", "", "1500\n", "", "", "", "", "2000\n", "2100\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(0+2200, 100), ylim = c(10, 80000), expand = FALSE, clip = "off") +
  labs(x = "Years (i.e., 2200 - year); logged", y = "Population (in millions); logged") +
  geom_hline(yintercept = 10) + geom_vline(xintercept = c(2200, 100)) +
  
  #Present phase
  annotate("segment", x =  1200, xend = 750, y = 6600, yend = 500, size = 0.25, colour = "#1D91C0", arrow = arrow(length = unit(0.2, "cm"))) +
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(2200), -log10(100))) + 
  annotate("label", colour = "#1D91C0", x = 1950, y = 10000, hjust = 0, vjust = 0.5, label = as.character(expression(paste(italic(P[T])," " == " ",over(paste(3, " ", 820) %*% 10^6,paste("ln",group("[",1.25 + italic(e)^paste((1980-italic(t))/25.5),"]")^0.716)),"  ,   ", t>+400))), parse = TRUE) +
  
  # Header
  annotate("rect", ymin = 80000, ymax = 25000, xmin = 100, xmax = 2200, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(2200, 100)))), y = exp(mean(log(c(80000, 25000)))), hjust = 0.5, vjust = 0.5, label = "Present phase", size = 13/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text.x = element_text(vjust = 0.5, hjust = 0.5),
        panel.grid.major.x = element_line(colour = "gray90"), panel.grid.minor.x = element_line(colour = "gray95"),
        panel.grid.major.y = element_line(colour = "gray90"), panel.grid.minor.y = element_blank(),
        plot.margin=unit(c(0.5, 0.5, 0.5, 0.5),"cm"))



# Saving the output
ggsave(filename = "04-figures and models/Fig1 - present phase.png", width = 17.5, height = 12.5, units = "cm", dpi = 600)





# __ FIGURE 2: Ancient period -----------------------------------------------

ggplot(subset(population, year < 400), 
       aes(x = year_min2200, y = baseline/1000000)) +
  
  geom_point(size = 1) +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.005, size = 0.2) +
  scale_x_continuous(breaks = c(1700, seq(2200, 12200, 1000), seq(22200, 102200, 10000), seq(102200, 1002200, 100000)),
                     labels = c("+500\n", "0(1)\n", "", "", "", "", "", "", "", "", "", "-10 000\n", "", "", "", "", "", "", "", "", "-100 000\n", 
                                "", "", "", "", "", "", "", "", "", "-1 000 000\n"),  
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(10^6+2200, 1700), ylim = c(0.001, 10000), expand = FALSE, clip = "off") +
  labs(x = "Years (i.e., 2200 - year); logged", y = "Population (in millions); logged") +
  geom_hline(yintercept = 0.001) + geom_vline(xintercept = c(1700, 10^6+2200)) +
  
  # Ancient phase
  annotate("segment", x = 8^6, xend = 7.5^6, y = 0.01, yend = 0.30, size = 0.25, colour = "#41AB5D", arrow = arrow(length = unit(0.2, "cm"))) +
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  annotate("label", colour = "#41AB5D", x = 9^6, y = 0.01, hjust = 0, vjust = 0.9, label = as.character(expression(paste(italic(P[T])," " == " ",over(paste(2, " ", 300) %*% 10^6,paste("ln",group("[","34,000" + italic(e)^paste((100-italic(t))/25.5),"]"))),"  ,   ", t<+400))), parse = TRUE) +
  
  # Header
  annotate("rect", xmin = 1002200, xmax = 1700, ymin = 1000, ymax = 10000, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1002200, 1700)))), y = exp(mean(log(c(1000, 10000)))), hjust = 0.5, vjust = 0.5, label = "Ancient phase", size = 13/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text.x = element_text(vjust = 0.5, hjust = 0.5),
        panel.grid.major.x = element_line(colour = "gray90"), panel.grid.minor.x = element_line(colour = "gray95"),
        panel.grid.major.y = element_line(colour = "gray90"), panel.grid.minor.y = element_blank(),
        plot.margin=unit(c(0.5, 0.5, 0.5, 0.5),"cm"))



# Saving the output
ggsave(filename = "04-figures and models/Fig2 - ancient phase.png", width = 17.5, height = 12.5, units = "cm", dpi = 600)





# __ FIGURE 3: Full figure ---------------------------------------------------------------

ggplot(population, aes(x = year_min2200, y = baseline/1000000)) +
  geom_vline(xintercept = c(200, 1200, 2200, 12200, 102200, 102200), color = "grey75") +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 
                                1700, 1800, 1900, 2000, 2100, 
                                seq(2200, 12200, 1000), seq(22200, 102200, 10000), seq(102200, 1002200, 100000)),
                     labels = c("", "2000\n", "", "", "", "", "", "", "", "", "", "1000\n", "", "", "", "", 
                                "", "", "", "", "", "0(1)\n", "", "", "", "", "", "", "", "", "", "-10 000\n", "", "", "", "", "", "", "", "", "-100 000\n", 
                                "", "", "", "", "", "", "", "", "", "-1 000 000\n"),  
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 30000), expand = FALSE, clip = "off") +
  labs(x = "Years (i.e., 2200 - year); logged", y = "Population (in millions); logged") +
  
  
  #Present phase
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  
  #Ancient phase
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text.x = element_text(vjust = 0.5, hjust = 0.5),
        panel.grid.major.x = element_line(colour = "gray90"), panel.grid.minor.x = element_line(colour = "gray95"),
        panel.grid.major.y = element_line(colour = "gray90"), panel.grid.minor.y = element_blank(),
        plot.margin=unit(c(0.5, 0.5, 0.5, 0.5),"cm"))

# Saving the output
ggsave(filename = "04-figures and models/Fig3 - full time range.png", width = 20, height = 12.5, units = "cm", dpi = 600)






# FIGURE 4: Average world population estimates ----------------------------

# Loading data from Table 1
table1 <- read_excel("01-data/Table1.xlsx")
table1[19, ] <- NA
table1$t <- as.numeric(table1$t)


# Defining functions
#E = exp((1980-x)/25.5)
f.PT <- function(x) 3.82/log(1.25+exp((1980-x)/25.5))^0.716
f.CT <- function(x) 10.94/log(2.65+exp((1980-x)/25.5))^0.852
f.TT <- function(x) (1.88/log(1+exp((1980-x)/25.5))^0.782)


ggplot(table1, aes(x = t, y = act.P)) +
  geom_hline(yintercept = c(0.1, 1, 10), color = "grey75") +
  stat_function(fun = f.CT, colour = "#41AB5D", xlim = c(400, 2100)) + annotate("text", x = 505, y = log10(2.58), hjust = 0, label = as.character(expression(paste(italic(C[T])))), parse = TRUE, color = "#41AB5D") +
  stat_function(fun = f.TT, colour = "#EF3B2C", xlim = c(400, 2100)) + annotate("text", x = 505, y = log10(1.25), hjust = 0, label = as.character(expression(paste(italic(T[T])))), parse = TRUE, color = "#EF3B2C") +
  stat_function(fun = f.PT, colour = "#1D91C0", xlim = c(400, 2100)) + annotate("text", x = 505, y = log10(1.78), hjust = 0, label = as.character(expression(paste(italic(P[T])))), parse = TRUE, color = "#1D91C0")  +
  geom_point() +
  scale_y_continuous(expression(paste("Population (", italic(P),"), Carrying capacity (", italic(C),"); in billion")),
                     breaks = c(0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15),
                     labels = c("0.01", "", "", "", "0.05", "", "", "", "", "0.1", "", "", "", "0.5", "", "", "", "", "1", "", "", "", "5", "", "", "", "", "10", "", "", "", "", "15"),
                     trans = "log10",
                     sec.axis = dup_axis(~./20, name = expression(paste("Technology (", italic(T),")")),)) +
  scale_x_continuous("Year") +
  coord_cartesian(xlim = c(0, 2100), ylim = c(0.05, 15), expand = FALSE, clip = "off") +
  theme(strip.background = element_rect(colour = "black"), 
        panel.background = element_blank(),
        axis.title.y.right = element_text(angle = 90),
        axis.ticks = element_blank(),
        panel.grid.major.x = element_line(colour = "gray90"), panel.grid.minor.x = element_line(colour = "gray95"),
        panel.grid.major.y = element_line(colour = "gray90"), panel.grid.minor.y = element_blank(),
        plot.margin=unit(c(0.5, 0.5, 0.5, 0.5),"cm"))

ggsave(filename = "04-figures and models/Fig4 - Average world population estimates.png", width = 17.5, height = 11.5, units = "cm", dpi = 600)


# Cleaning up
rm(f.PT, f.CT, f.TT)





# FIGURE A1: Individual figures for each dataset --------------------------


# __ Generating individual figures for each dataset -----------------------

# Taagepera (1979)
Taagepera_1979 <-
  ggplot(subset(population, source == "Taagepera (1979)"), 
         aes(x = year_min2200, y = baseline/1000000)) +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(10^6+2200, 10^5+2200, 10^4+2200, 4200, 2200, 1200, 700, 200, 100, 50),
                     labels = c("-1 000 000\n", "-100 000\n", "-10 000\n", "\n-2000", "0(1)\n", "\n1000", "1500\n", "2000\n", "2100\n", "2150\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 100000), expand = FALSE, clip = "off") +
  labs(x = NULL, y = NULL) +
  geom_vline(xintercept = c(10^6, 100)) +
  geom_hline(yintercept = 0.001) +
  
  # Functions
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  # Header
  annotate("rect", ymin = 100000, ymax = 14000, xmin = 100, xmax = 10^6, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1000000, 100)))), y = exp(mean(log(c(100000, 14000)))), hjust = 0.5, vjust = 0.5, label = "Taagepera (1979)", size = 11/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text.y = element_text(vjust = 0.5, hjust = 0.5),
        axis.text.x = element_blank(),
        panel.grid.major.x = element_line(colour = "gray90"), 
        panel.grid.major.y = element_line(colour = "gray90"), 
        plot.margin=unit(c(0.05, 0.35, 0.05, 0.35),"cm"))



# HYDE 3.2
hyde_3_2 <-
  ggplot(subset(population, source == "HYDE 3.2"), 
         aes(x = year_min2200, y = baseline/1000000)) +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(10^6+2200, 10^5+2200, 10^4+2200, 4200, 2200, 1200, 700, 200, 100, 50),
                     labels = c("-1 000 000\n", "-100 000\n", "-10 000\n", "\n-2000", "0(1)\n", "\n1000", "1500\n", "2000\n", "2100\n", "2150\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 100000), expand = FALSE, clip = "off") +
  labs(x = NULL, y = NULL) +
  geom_vline(xintercept = c(10^6, 100)) +
  geom_hline(yintercept = 0.001) +
  
  # Functions
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  # Header
  annotate("rect", ymin = 100000, ymax = 14000, xmin = 100, xmax = 10^6, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1000000, 100)))), y = exp(mean(log(c(100000, 14000)))), hjust = 0.5, vjust = 0.5, label = "HYDE 3.2 (History database of the Global Environment)", size = 11/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        panel.grid.major.x = element_line(colour = "gray90"), 
        panel.grid.major.y = element_line(colour = "gray90"), 
        plot.margin=unit(c(0.05, 0.35, 0.05, 0.35),"cm"))



# McEvedy & Jones (1978)
McEvedy_Jones_1978 <-
  ggplot(subset(population, source == "McEvedy & Jones (1978)"), 
         aes(x = year_min2200, y = baseline/1000000)) +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(10^6+2200, 10^5+2200, 10^4+2200, 4200, 2200, 1200, 700, 200, 100, 50),
                     labels = c("-1 000 000\n", "-100 000\n", "-10 000\n", "\n-2000", "0(1)\n", "\n1000", "1500\n", "2000\n", "2100\n", "2150\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 100000), expand = FALSE, clip = "off") +
  labs(x = NULL, y = NULL) +
  geom_vline(xintercept = c(10^6, 100)) +
  geom_hline(yintercept = 0.001) +
  
  # Functions
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  # Header
  annotate("rect", ymin = 100000, ymax = 14000, xmin = 100, xmax = 10^6, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1000000, 100)))), y = exp(mean(log(c(100000, 14000)))), hjust = 0.5, vjust = 0.5, label = "McEvedy & Jones (1978)", size = 11/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text.x = element_blank(),
        panel.grid.major.x = element_line(colour = "gray90"), 
        panel.grid.major.y = element_line(colour = "gray90"), 
        plot.margin=unit(c(0.05, 0.35, 0.05, 0.35),"cm"))



# Thomlinson (1975)
Thomlinson_1975 <-
  ggplot(subset(population, source == "Thomlinson (1975)"), 
         aes(x = year_min2200, y = baseline/1000000)) +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(10^6+2200, 10^5+2200, 10^4+2200, 4200, 2200, 1200, 700, 200, 100, 50),
                     labels = c("-1 000 000\n", "-100 000\n", "-10 000\n", "\n-2000", "0(1)\n", "\n1000", "1500\n", "2000\n", "2100\n", "2150\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 100000), expand = FALSE, clip = "off") +
  labs(x = NULL, y = NULL) +
  geom_vline(xintercept = c(10^6, 100)) +
  geom_hline(yintercept = 0.001) +
  
  # Functions
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  # Header
  annotate("rect", ymin = 100000, ymax = 14000, xmin = 100, xmax = 10^6, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1000000, 100)))), y = exp(mean(log(c(100000, 14000)))), hjust = 0.5, vjust = 0.5, label = "Thomlinson (1975)", size = 11/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        panel.grid.major.x = element_line(colour = "gray90"), 
        panel.grid.major.y = element_line(colour = "gray90"), 
        plot.margin=unit(c(0.05, 0.35, 0.05, 0.35),"cm"))



# Population Reference Bureau (PRB)
PRB <-
  ggplot(subset(population, source == "Population Reference Bureau (PRB)"), 
         aes(x = year_min2200, y = baseline/1000000)) +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(10^6+2200, 10^5+2200, 10^4+2200, 4200, 2200, 1200, 700, 200, 100, 50),
                     labels = c("-1 000 000\n", "-100 000\n", "-10 000\n", "\n-2000", "0(1)\n", "\n1000", "1500\n", "2000\n", "2100\n", "2150\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 100000), expand = FALSE, clip = "off") +
  labs(x = NULL, y = NULL) +
  geom_vline(xintercept = c(10^6, 100)) +
  geom_hline(yintercept = 0.001) +
  
  # Functions
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  # Header
  annotate("rect", ymin = 100000, ymax = 14000, xmin = 100, xmax = 10^6, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1000000, 100)))), y = exp(mean(log(c(100000, 14000)))), hjust = 0.5, vjust = 0.5, label = "Population Reference Bureau (PRB)", size = 11/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text.x = element_blank(),
        panel.grid.major.x = element_line(colour = "gray90"), 
        panel.grid.major.y = element_line(colour = "gray90"), 
        plot.margin=unit(c(0.05, 0.35, 0.05, 0.35),"cm"))



# Durand (1974)
Durand_1974 <-
  ggplot(subset(population, source == "Durand (1974)"), 
         aes(x = year_min2200, y = baseline/1000000)) +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(10^6+2200, 10^5+2200, 10^4+2200, 4200, 2200, 1200, 700, 200, 100, 50),
                     labels = c("-1 000 000\n", "-100 000\n", "-10 000\n", "\n-2000", "0(1)\n", "\n1000", "1500\n", "2000\n", "2100\n", "2150\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 100000), expand = FALSE, clip = "off") +
  labs(x = NULL, y = NULL) +
  geom_vline(xintercept = c(10^6, 100)) +
  geom_hline(yintercept = 0.001) +
  
  # Functions
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  # Header
  annotate("rect", ymin = 100000, ymax = 14000, xmin = 100, xmax = 10^6, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1000000, 100)))), y = exp(mean(log(c(100000, 14000)))), hjust = 0.5, vjust = 0.5, label = "Durand (1974)", size = 11/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        panel.grid.major.x = element_line(colour = "gray90"), 
        panel.grid.major.y = element_line(colour = "gray90"), 
        plot.margin=unit(c(0.05, 0.35, 0.05, 0.35),"cm"))



# United Nations
UN <-
  ggplot(subset(population, source == "United Nations"), 
         aes(x = year_min2200, y = baseline/1000000)) +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(10^6+2200, 10^5+2200, 10^4+2200, 4200, 2200, 1200, 700, 200, 100, 50),
                     labels = c("-1 000 000\n", "-100 000\n", "-10 000\n", "\n-2000", "0(1)\n", "\n1000", "1500\n", "2000\n", "2100\n", "2150\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 100000), expand = FALSE, clip = "off") +
  labs(x = NULL, y = NULL) +
  geom_vline(xintercept = c(10^6, 100)) +
  geom_hline(yintercept = 0.001) +
  
  # Functions
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  # Header
  annotate("rect", ymin = 100000, ymax = 14000, xmin = 100, xmax = 10^6, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1000000, 100)))), y = exp(mean(log(c(100000, 14000)))), hjust = 0.5, vjust = 0.5, label = "United Nations Department of Economic and Social Affairs", size = 11/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text.x = element_text(vjust = 0.5, hjust = 0.5),
        panel.grid.major.x = element_line(colour = "gray90"), 
        panel.grid.major.y = element_line(colour = "gray90"), 
        plot.margin=unit(c(0.05, 0.35, 0.05, 0.35),"cm"))



# Maddison (2008)
Maddison_2008 <-
  ggplot(subset(population, source == "Maddison (2008)"), 
         aes(x = year_min2200, y = baseline/1000000)) +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(10^6+2200, 10^5+2200, 10^4+2200, 4200, 2200, 1200, 700, 200, 100, 50),
                     labels = c("-1 000 000\n", "-100 000\n", "-10 000\n", "\n-2000", "0(1)\n", "\n1000", "1500\n", "2000\n", "2100\n", "2150\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 100000), expand = FALSE, clip = "off") +
  labs(x = NULL, y = NULL) +
  geom_vline(xintercept = c(10^6, 100)) +
  geom_hline(yintercept = 0.001) +
  
  # Functions
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  # Header
  annotate("rect", ymin = 100000, ymax = 14000, xmin = 100, xmax = 10^6, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1000000, 100)))), y = exp(mean(log(c(100000, 14000)))), hjust = 0.5, vjust = 0.5, label = "Maddison (2008)", size = 11/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text.y = element_blank(),
        axis.text.x = element_text(vjust = 0.5, hjust = 0.5),
        panel.grid.major.x = element_line(colour = "gray90"), 
        panel.grid.major.y = element_line(colour = "gray90"), 
        plot.margin=unit(c(0.05, 0.2, 0.05, 0),"cm"))




# SECOND SET OF FIGURES
# Tanton (1994)
Tanton_1994 <-
  ggplot(subset(population, source == "Tanton (1994)"), 
         aes(x = year_min2200, y = baseline/1000000)) +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(10^6+2200, 10^5+2200, 10^4+2200, 4200, 2200, 1200, 700, 200, 100, 50),
                     labels = c("-1 000 000\n", "-100 000\n", "-10 000\n", "\n-2000", "0(1)\n", "\n1000", "1500\n", "2000\n", "2100\n", "2150\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 100000), expand = FALSE, clip = "off") +
  labs(x = NULL, y = NULL) +
  geom_vline(xintercept = c(10^6, 100)) +
  geom_hline(yintercept = 0.001) +
  
  # Functions
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  # Header
  annotate("rect", ymin = 100000, ymax = 14000, xmin = 100, xmax = 10^6, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1000000, 100)))), y = exp(mean(log(c(100000, 14000)))), hjust = 0.5, vjust = 0.5, label = "Tanton (1994)", size = 11/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text.x = element_blank(),
        panel.grid.major.x = element_line(colour = "gray90"), 
        panel.grid.major.y = element_line(colour = "gray90"), 
        plot.margin=unit(c(0.05, 0.35, 0.05, 0.35),"cm"))



# Biraben (1980)
Biraben_1980 <-
  ggplot(subset(population, source == "Biraben (1980)"), 
         aes(x = year_min2200, y = baseline/1000000)) +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(10^6+2200, 10^5+2200, 10^4+2200, 4200, 2200, 1200, 700, 200, 100, 50),
                     labels = c("-1 000 000\n", "-100 000\n", "-10 000\n", "\n-2000", "0(1)\n", "\n1000", "1500\n", "2000\n", "2100\n", "2150\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 100000), expand = FALSE, clip = "off") +
  labs(x = NULL, y = NULL) +
  geom_vline(xintercept = c(10^6, 100)) +
  geom_hline(yintercept = 0.001) +
  
  # Functions
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  # Header
  annotate("rect", ymin = 100000, ymax = 14000, xmin = 100, xmax = 10^6, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1000000, 100)))), y = exp(mean(log(c(100000, 14000)))), hjust = 0.5, vjust = 0.5, label = "Biraben (1980)", size = 11/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        panel.grid.major.x = element_line(colour = "gray90"), 
        panel.grid.major.y = element_line(colour = "gray90"), 
        plot.margin=unit(c(0.05, 0.35, 0.05, 0.35),"cm"))



# Clark (1967)
Clark_1967 <-
  ggplot(subset(population, source == "Clark (1967)"), 
         aes(x = year_min2200, y = baseline/1000000)) +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(10^6+2200, 10^5+2200, 10^4+2200, 4200, 2200, 1200, 700, 200, 100, 50),
                     labels = c("-1 000 000\n", "-100 000\n", "-10 000\n", "\n-2000", "0(1)\n", "\n1000", "1500\n", "2000\n", "2100\n", "2150\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 100000), expand = FALSE, clip = "off") +
  labs(x = NULL, y = NULL) +
  geom_vline(xintercept = c(10^6, 100)) +
  geom_hline(yintercept = 0.001) +
  
  # Functions
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  # Header
  annotate("rect", ymin = 100000, ymax = 14000, xmin = 100, xmax = 10^6, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1000000, 100)))), y = exp(mean(log(c(100000, 14000)))), hjust = 0.5, vjust = 0.5, label = "Clark (1967)", size = 11/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text.x = element_blank(),
        panel.grid.major.x = element_line(colour = "gray90"), 
        panel.grid.major.y = element_line(colour = "gray90"), 
        plot.margin=unit(c(0.05, 0.35, 0.05, 0.35),"cm"))



# Gapminder 6.0
gapminder_v6 <-
  ggplot(subset(population, source == "Gapminder 6.0"), 
         aes(x = year_min2200, y = baseline/1000000)) +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(10^6+2200, 10^5+2200, 10^4+2200, 4200, 2200, 1200, 700, 200, 100, 50),
                     labels = c("-1 000 000\n", "-100 000\n", "-10 000\n", "\n-2000", "0(1)\n", "\n1000", "1500\n", "2000\n", "2100\n", "2150\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 100000), expand = FALSE, clip = "off") +
  labs(x = NULL, y = NULL) +
  geom_vline(xintercept = c(10^6, 100)) +
  geom_hline(yintercept = 0.001) +
  
  # Functions
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  # Header
  annotate("rect", ymin = 100000, ymax = 14000, xmin = 100, xmax = 10^6, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1000000, 100)))), y = exp(mean(log(c(100000, 14000)))), hjust = 0.5, vjust = 0.5, label = "Gapminder (Version 6.0)", size = 11/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        panel.grid.major.x = element_line(colour = "gray90"), 
        panel.grid.major.y = element_line(colour = "gray90"), 
        plot.margin=unit(c(0.05, 0.35, 0.05, 0.35),"cm"))



# UN Population Prospects (2019)
UN_PP_2019 <-
  ggplot(subset(population, source == "UN Population Prospects (2019)"), 
         aes(x = year_min2200, y = baseline/1000000)) +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(10^6+2200, 10^5+2200, 10^4+2200, 4200, 2200, 1200, 700, 200, 100, 50),
                     labels = c("-1 000 000\n", "-100 000\n", "-10 000\n", "\n-2000", "0(1)\n", "\n1000", "1500\n", "2000\n", "2100\n", "2150\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 100000), expand = FALSE, clip = "off") +
  labs(x = NULL, y = NULL) +
  geom_vline(xintercept = c(10^6, 100)) +
  geom_hline(yintercept = 0.001) +
  
  # Functions
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  # Header
  annotate("rect", ymin = 100000, ymax = 14000, xmin = 100, xmax = 10^6, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1000000, 100)))), y = exp(mean(log(c(100000, 14000)))), hjust = 0.5, vjust = 0.5, label = "UN Population Prospects (2019)", size = 11/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text.x = element_text(vjust = 0.5, hjust = 0.5),
        panel.grid.major.x = element_line(colour = "gray90"), 
        panel.grid.major.y = element_line(colour = "gray90"), 
        plot.margin=unit(c(0.05, 0.35, 0.05, 0.35),"cm"))



# US Census Bureau (2017)
US_CB_2017 <-
  ggplot(subset(population, source == "US Census Bureau (2017)"), 
         aes(x = year_min2200, y = baseline/1000000)) +
  geom_point() +
  geom_errorbar(aes(x = year_min2200,  ymin = lower/1000000,  ymax = upper/1000000), width = 0.05) +
  scale_x_continuous(breaks = c(10^6+2200, 10^5+2200, 10^4+2200, 4200, 2200, 1200, 700, 200, 100, 50),
                     labels = c("-1 000 000\n", "-100 000\n", "-10 000\n", "\n-2000", "0(1)\n", "\n1000", "1500\n", "2000\n", "2100\n", "2150\n"),
                     trans = reverselog_trans(10), minor_breaks = FALSE) +
  scale_y_continuous(breaks = c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000), 
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100", "1 000", "10 000"), 
                     trans = "log10", minor_breaks = FALSE) +
  coord_cartesian(xlim = c(1002200, 100), ylim = c(0.001, 100000), expand = FALSE, clip = "off") +
  labs(x = NULL, y = NULL) +
  geom_vline(xintercept = c(10^6, 100)) +
  geom_hline(yintercept = 0.001) +
  
  # Functions
  stat_function(fun = f.PT.new.1, colour = "#1D91C0", xlim = c(-log10(1800), -log10(100))) + 
  stat_function(fun = f.PT.old.1, colour = "#41AB5D", xlim = c(-4.30103, -3.230449)) +
  stat_function(fun = f.PT.old.2, colour = "#41AB5D", xlim = c(-6, -4.30103)) +
  
  # Header
  annotate("rect", ymin = 100000, ymax = 14000, xmin = 100, xmax = 10^6, colour = "black", fill = "grey90") +
  annotate("text", x = exp(mean(log(c(1000000, 100)))), y = exp(mean(log(c(100000, 14000)))), hjust = 0.5, vjust = 0.5, label = "US Census Bureau (2017)", size = 11/.pt) +
  
  theme(panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.text.y = element_blank(),
        axis.text.x = element_text(vjust = 0.5, hjust = 0.5),
        panel.grid.major.x = element_line(colour = "gray90"), 
        panel.grid.major.y = element_line(colour = "gray90"), 
        plot.margin=unit(c(0.05, 0.35, 0.05, 0.35),"cm"))







# __ Figure A1 - PART 1 ---------------------------------------------------
annotate_figure(
  ggarrange(Taagepera_1979, hyde_3_2, 
            McEvedy_Jones_1978, Thomlinson_1975, 
            PRB, Durand_1974, 
            UN, Maddison_2008,
            nrow = 4, ncol = 2,
            align = "hv"),
  left = textGrob(expression(paste("Population (in millions); logged")), rot = 90, hjust = 0.5, vjust = 0.5, gp = gpar(fontsize = 11, col = "black", fill = "white")),
  bottom = textGrob(expression(paste("Years (i.e., 2200 - year); logged")), hjust = 0.5, vjust = 0.5, gp = gpar(fontsize = 11, col = "black", fill = "white")),
  right = ""
) + theme(plot.background = element_rect(fill = "white", color = "transparent"))




# Saving the output
ggsave(filename = "04-figures and models/FigA1a - individual data sources.png", width = 26.5, height = 25, units = "cm", dpi = 600)




# __ Figure A1 - PART 2 ---------------------------------------------------
annotate_figure(
  ggarrange(Tanton_1994, Biraben_1980,
            Clark_1967, gapminder_v6, 
            UN_PP_2019, US_CB_2017, 
            nrow = 3, ncol = 2, 
            align = "hv"),
  left = textGrob(expression(paste("Population (in millions); logged")), rot = 90, hjust = 0.5, vjust = 0.5, gp = gpar(fontsize = 11, col = "black", fill = "white")),
  bottom = textGrob(expression(paste("Years (i.e., 2200 - year); logged")), hjust = 0.5, vjust = 0.5, gp = gpar(fontsize = 11, col = "black", fill = "white")),
  right = ""
) + theme(plot.background = element_rect(fill = "white", color = "transparent"))




# Saving the output
ggsave(filename = "04-figures and models/FigA1b - individual data sources.png", width = 26.5, height = 18.75, units = "cm", dpi = 600)



# Cleaning up
rm(hyde_3_2, gapminder_v6, UN_PP_2019, McEvedy_Jones_1978, Thomlinson_1975, PRB, Durand_1974, 
   UN, Maddison_2008, Tanton_1994, Biraben_1980, Clark_1967, US_CB_2017, Taagepera_1979)

rm(f.PT.new.1, f.PT.new.2, f.PT.old.1, f.PT.old.2)
rm(reverselog_trans)



# FIGURE A2: Redundancy factor f[PT/U] vs PT/U ----------------------------
table1$PTU <- table1$PT/table1$U
table1$fPTU <- table1$"f[PT/U]"

f.PT <- function(x) (1+exp(75*(x-0.13)))^-0.023

ggplot(table1, aes(x = PTU, y = fPTU, label = t)) +
  geom_point() +
  stat_function(fun = f.PT, size = 1, colour = "#1D91C0") +
  scale_y_continuous(breaks = c(0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1), 
                     labels = c("0.2", "", "", "0.5", "", "", "", "", "1"), trans = "log10") +
  scale_x_continuous(breaks = c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1), 
                     labels = c("0", "", "", "", "", "0.5", "", "", "", "", "1")) +
  coord_cartesian(xlim = c(0, 1), ylim = c(0.2, 1.015), expand = FALSE, clip = "off") +
  geom_segment(aes(x = 0, xend = 0, y = 0.2, yend = 1)) +
  geom_segment(aes(x = 0, xend = 1, y = 0.2, yend = 0.2)) +
  labs(x = expression(paste(italic(P[T])/italic(U))), y = expression(paste(italic(f),"[",italic(P[T]/U),"]; logged"))) +
  annotate("label", x = 0.98, y = 0.98, hjust = 1, vjust = 1,
           label = as.character(expression(paste(italic(f),"[",italic(P[T]/U),"]" == "[",1 + italic(e)^75(P[T]/U-0.13),"]"^-0.023))), 
           parse = TRUE) +
  theme(legend.position = "none", panel.border = element_blank(),
        strip.background = element_blank(),
        panel.grid.major.x = element_line(colour = "gray95"),
        panel.grid.minor.x = element_line(colour = "gray95"),
        panel.grid.major.y = element_line(colour = "gray95"),
        panel.grid.minor.y = element_line(colour = "gray95"),
        panel.background = element_blank(),
        plot.margin=unit(c(0.5, 0.5, 0.5, 0.5),"cm")) 
  
# Saving the output
ggsave(filename = "04-figures and models/FigA2 - Redundancy factor.png", width = 17.5, height = 12.5, units = "cm", dpi = 600)


# Final cleanup
rm(table1, f.PT)
rm(population)
